﻿using System;
using System.Collections.Generic;
using System.Text;
using VerticalProcesses;

namespace VerticalProcesses
{
    /// <summary>
    /// Basic interface to calculate engergy input onto the snowpack
    /// </summary>
    public abstract class SnowpackEnergyInputEq
    {
        //variables
        protected double TimeStepHrs = 1.0;
        
        //properties
        public double SnowpackEnergyInput { get; set; }              
        public abstract double CalcSnowpackEnergyInput(DateTime TimeStep, double gruTa, double gruP, double gruShortRad, double tempSurf, double crownClosure);


    }
}
# region Implementions of different equations to calculate engergy input onto the snowpack
/// <summary>
/// Add your snowpack implemention here
/// </summary>


public class JostMoore : SnowpackEnergyInputEq
{
    //global parameters set in constructor
    private double AF;               // temperature component of snowmelt factor at clearcut
    private double BF;               // temperature component of snowmelt factor at crown closure = 1
    private double CF;               // radiation component of snowmelt factor
    private double T0;               // threshold temperature for snowmelt
    private double CCRF;             // reduction factor for cold content   

    
    /// <summary>
    /// Calculates snowpack energy input according to a 
    /// modification of the distriuted temperature index
    /// formulation propose by Hock (1999)
    /// </summary>
    /// <param name="TimeStepHrs [h]"></param>
    /// <param name="thrfRain [mm]"></param>
    /// <param name="Radiation [W/m2]"></param>
    /// <param name="gruTa[ dC]"></param>
    /// <returns>Engergy input to snowpack  [MJ/timestep/m2]</returns>
    public override double CalcSnowpackEnergyInput(DateTime TimeStep, double gruTa, double gruRain, double radiation, double tempSurf, double crownClosure)
    {

        //provide an option to calculate cold content
        double MeltMJ;
        //radiation = radiation * 3600.0 * TimeStepHrs / 1.0E6; //convert W to MJ (60 [sec/min] * 60 [min/h])                      
        double mf = AF * (1 - crownClosure) + BF * crownClosure;
        //double mf = AF * TimeStepHrs;
        //double mf = AF * Math.Exp(-BF * crownClosure);
        if (gruTa < 0.0)
        {
            if (gruTa < tempSurf)
            {

                SnowpackEnergyInput = CCRF * (gruTa - tempSurf);
                //SnowpackEnergyInput = CCRF * (mf - CF * radiation) * (gruTa - tempSurf);
                //this.SAFCC = 1 + ALPHAC * Math.Sin(slope * Math.PI / 180.0) * Math.Cos(aspect * Math.PI / 180.0);

            }
            
                else 
            {
                SnowpackEnergyInput = (mf + CF * radiation) * (gruTa - tempSurf);
                //SnowpackEnergyInput = Math.Max(mf, CF * radiation) * (gruTa - tempSurf);
            }
        }
        else
        {
            double rainEnergyInput = gruRain * Math.Max(gruTa, 0.0) * Constants.CvLq;
             MeltMJ = (mf + CF * radiation) * (gruTa - T0);
            //MeltMJ = Math.Max(mf, CF * radiation) * (gruTa - T0);
            SnowpackEnergyInput = MeltMJ + rainEnergyInput;
        }
        
        return SnowpackEnergyInput;
    }
    public JostMoore(double AF, double BF, double CF, double CCRF, double T0)
    {
        this.AF = AF;
        this.BF = BF;
        this.CF = CF;
        this.CCRF = CCRF;
        this.T0 = T0;
    }


}
public class DegreeDay : SnowpackEnergyInputEq
{
    private double AF;               // temperature component of snowmelt factor at clearcut
    private double BF;               // temperature component of snowmelt factor at crown closure = 1
    private double T0;               // threshold temperature for snowmelt
    private double CCRF;             // reduction factor for cold content 
    
    /// <summary>
    /// Simple degree day model
    /// </summary>
    /// <param name="gruTa"></param>
    /// <param name="gruRain"></param>
    /// <param name="tempSurf"></param>
    /// <param name="crownClosure"></param>
    /// <returns>Engergy input to snowpack  [MJ/timestep/m2]</returns>
    public override double CalcSnowpackEnergyInput(DateTime TimeStep, double gruTa, double gruRain, double radiation, double tempSurf, double crownClosure)
    {

        //provide an option to calculate cold content
        double MeltMJ;
        double mf = (AF * (1 - crownClosure) + BF * crownClosure) * TimeStepHrs;


        if (gruTa < 0.0)
        {
            if (gruTa < tempSurf)
            {
                SnowpackEnergyInput = CCRF * (gruTa - tempSurf);
            }
            else
            {
                SnowpackEnergyInput = mf * (gruTa - tempSurf);
            }
        }
        else
        {
            double rainEnergyInput = gruRain * Math.Max(gruTa, 0.0) * Constants.CvLq;
            MeltMJ = mf * (gruTa - T0);
            SnowpackEnergyInput = MeltMJ + rainEnergyInput;
        }

        return SnowpackEnergyInput;
    }
    public DegreeDay(double AF, double BF, double CF, double CCRF, double T0)
    {
        this.AF = AF;
        this.BF = BF;
        this.CCRF = CCRF;
        this.T0 = T0;
    }

}
public class Kustas : SnowpackEnergyInputEq
{
    //global parameters set in constructor
    private double AF;               // temperature component of snowmelt factor at clearcut
    private double BF;               // temperature component of snowmelt factor at crown closure = 1
    private double CF;               // radiation component of snowmelt factor
    private double T0;               // threshold temperature for snowmelt
    private double CCRF;             // reduction factor for cold content   


    /// <summary>
    /// Calculates snowpack energy input according to Kustas (1994)
    /// </summary>
    /// <param name="TimeStepHrs [h]"></param>
    /// <param name="thrfRain [mm]"></param>
    /// <param name="Radiation [W/m2]"></param>
    /// <param name="gruTa[ dC]"></param>
    /// <returns>Engergy input to snowpack  [MJ/timestep/m2]</returns>
    public override double CalcSnowpackEnergyInput(DateTime TimeStep, double gruTa, double gruRain, double radiation, double tempSurf, double crownClosure)
    {

        //provide an option to calculate cold content
        double MeltMJ;
        //radiation = radiation * 3600.0 * TimeStepHrs / 1.0E6; //convert W to MJ (60 [sec/min] * 60 [min/h])                               
        double mf = (AF * (1 - crownClosure) + BF * crownClosure) * TimeStepHrs;
       

        if (gruTa < 0.0)
        {
            if (gruTa < tempSurf)
            {
                SnowpackEnergyInput = CCRF * (gruTa - tempSurf);
            }
            else
            {
                SnowpackEnergyInput = mf * (gruTa - tempSurf) + CF * radiation;
            }
        }
        else
        {
            double rainEnergyInput = gruRain * Math.Max(gruTa, 0.0) * Constants.CvLq;
            MeltMJ = mf * (gruTa - T0) + CF * radiation;
            SnowpackEnergyInput = MeltMJ + rainEnergyInput;
        }

        return SnowpackEnergyInput;
    }
    public Kustas(double AF, double BF, double CF, double CCRF, double T0)
    {
        this.AF = AF;
        this.BF = BF;
        this.CF = CF;
        this.CCRF = CCRF;
        this.T0 = T0;
    }


}
public class HBVEC : SnowpackEnergyInputEq
{
    //global parameters set in constructor
    private double AF;               // temperature component of snowmelt factor at clearcut
    private double BF;               // temperature component of snowmelt factor at crown closure = 1
    private double CF;               // radiation component of snowmelt factor
    private double T0;               // threshold temperature for snowmelt
    private double CCRF;             // reduction factor for cold content   
    private double SAF;              // slope aspect adjustment factor
    private double SAFCC;            // slope aspect adjustment factor for cold content

    /// <summary>
    /// Calculates snowpack energy input according to Kustas (1994)
    /// </summary>
    /// <param name="TimeStepHrs [h]"></param>
    /// <param name="thrfRain [mm]"></param>
    /// <param name="Radiation [W/m2]"></param>
    /// <param name="gruTa[ dC]"></param>
    /// <returns>Engergy input to snowpack  [MJ/timestep/m2]</returns>
    public override double CalcSnowpackEnergyInput(DateTime TimeStep, double gruTa, double gruRain, double radiation, double tempSurf, double crownClosure)
    {

        //provide an option to calculate cold content
        double MeltMJ;
        //variation of basemelt factor between solstices
        double daysInAYear = 365.25; // Average number of days in a year
        DateTime lastWintSolst = new DateTime(TimeStep.Year - 1, 12, 21, 12, 0, 0);
        TimeSpan daysSinceLastWintSolst = TimeStep.Subtract(lastWintSolst);//((TimeSpan)(dt.Subtract(winterSolstice))).TotalHours;
        double mf = (AF + 0.5 * BF * (1 - Math.Cos(2 * Math.PI * daysSinceLastWintSolst.Days / daysInAYear))) * TimeStepHrs;
        double mrf = 1 - crownClosure + crownClosure * CF;        

        if (gruTa < 0.0)
        {
            if (gruTa < tempSurf)
            {
                SnowpackEnergyInput = CCRF * (gruTa - tempSurf);
            }
            else
            {
                SnowpackEnergyInput = mf * mrf * SAF * (gruTa - tempSurf);
            }

        }
        else
        {
            double rainEnergyInput = gruRain * Math.Max(gruTa, 0.0) * Constants.CvLq;
            MeltMJ = mf * mrf * SAF * (gruTa - T0);
            SnowpackEnergyInput = MeltMJ + rainEnergyInput;
        }
        return SnowpackEnergyInput;
    }
    public HBVEC(double AF, double BF, double CF, double CCRF, double T0, double ALPHAC, double slope, double aspect)
    {
        this.AF = AF;
        this.BF = BF;
        this.CF = CF;
        this.CCRF = CCRF;
        this.T0 = T0;
        this.SAF = 1 - ALPHAC * Math.Sin(slope * Math.PI / 180.0) * Math.Cos(aspect * Math.PI / 180.0);
        this.SAFCC = 1 + ALPHAC * Math.Sin(slope * Math.PI / 180.0) * Math.Cos(aspect * Math.PI / 180.0);
    }


}
#endregion